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Abstract 

The recent discovery of the Higgs boson at the LHC attracts mnch attention to the 
precise calcnlation of its prodnction cross section in qnantnm chromodynamics. In 
this work, we discuss the calculation of soft triple-emission phase space integral, 
which is an essential ingredient in the recently calculated soft-virtual corrections 
to Higgs boson production at next-to-next-to-next-to-leading order. The main 
techniques used this calculation are method of differential equation for Feynman 
integral, and integration of harmonic polylogarithms. 


1 Introduction 


To maximize the physics outcome of the LHC program, theoretical predictions for signal and 
background processes have to match the unprecedented accuracy of experimental measurements. 
For QCD initiated processes, a standard approach to improve the accuracy of theoretical predic¬ 
tion is by calculating perturbative radiative corrections beyond Leading Order (LO) approxima¬ 
tion. At Next-to-Leading Order (NLO), such an calculation consists of virtual correction, which 
involves virtual parton circulated in the loop, and real correction, which involves real on-shell 
parton emission. When considering corrections beyond NLO, mixed real-virtual correction is 
also required. In the past few decades, signihcant efforts have been made to develop analytic 
techniques for loop amplitude calculation. In particular, a variety of methods have been devel¬ 
oped for the analytic calculation of loop integral, cf. [T] and references therein. For phase space 
integral, techniques for analytic calculation are considerably less developed. An obvious reason 
is that experimental measurements often involve complicated phase space cuts, which are very 
difficult to implement analytically. The standard approach is to perform the phase space integral 
numerically, keeping in mind that appropriate subtraction terms have to be added in order to 
render the integral hnite Eli. In that case analytic inclusive phase space integrals can be used 
in the construction of infrared subtraction terms for Next-to-Next-to-Leading Order (NNLO) 
calculation |1]. Moreover, when appropriate subtraction method was not yet available, analytic 
integral was the only method for obtaining cross section or distribution at NNLO at hadron 
collider [3 El [3 [HI i. 

Recently, an important step towards the calculation of Higgs production at Next-to-Next- 
to-Next-to-Leading Order (N^LO) has been taken in Refs. mm, in which the soft-virtual 
corrections are obtained. The soft-virtual corrections are cross section in the threshold limit, 
where radiation in the hnal state is constrained to have low energy. Since the radiation energy is 
low, calculation can be done in the eikonal limit, where QCD partons are emitted from light-like 
Wilson lines which parameterize the directions of incoming partons. Calculation of soft-virtual 
corrections at N^LO requires the calculation of pure virtual corrections through to three loops, 
and soft phase space integrals with one, two, and three partons in the hnal state. For soft phase 
space integrals, the ones with three partons in the hnal state are the most challenging one. Such 
triple soft phase space integrals are hrst calculated in the pioneering work [12], using various 
powerful techiques for Feynman integrals, ranging from Melin-Barnes transformation to symbols 
and coproduct. 

In Ref. mi, a completely diherent approach is taken. Auxiliary integrals are introduced in the 
intermediate steps, as well as nontrivial dependence on internal auxiliary scales. The auxiliary 
integrals are then solved by the method of diherential equation for Feynman integral. Finally the 
auxiliary scales are integrated over and the desired triple soft phase space integrals are obtained. 
This approach makes heavy use of the property of scale invariance for eikonal integral. It is 
worth mentioning that the triple soft phase space integrals are single scale integrals with trivial 
dependence on the single scale, therefore prohibit direct application of the method of differential 
equation. By introducing the auxiliary integrals, additional nontrivial scales are introduced into 
the problem, and nontrivial system of differential equation can be derived. While differential 
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equation approach to single scale integrals have been employed before in different context [laiii] 
by introducing external scales, we hnd that auxiliary scales being introduced internally in the 
current calculation is particularly interesting, and significantly simplifies the calculation at all 
stages. In this paper we document the detail of such a calculation, in the hope that it will be 
useful in other problems. 

This paper is organized as follows. In Sec. [21 we explain the method in detail, using the 
NNLO double soft emission phase space integral as example. In Sec. El we apply this method 
to triple soft phase space integrals. In Sec. 0] we give a brief conclusion and outlook. Detail for 
solving the system of differential equation for triple soft integrals is given in Appendix. [^ 

2 Differential equation approach to single scale soft phase 
space integral 

2.1 Introducing the method 

The method of differential equation [T^ [T6l [T7] is powerful in the calculation of Feynman inte¬ 
grals with multiple scales. For single scale problem, the dependence on the scale is trivial from 
dimensional analysis. The nontrivial part of the integral is a function of dimensional regulariza¬ 
tion parameter only 0. Thus the method of differential equation can not be applied directly. In 
some cases, one can introduce external auxiliary scales such that nontrivial differential equation 
for Feynman integral can be derived. The integrals with additional external auxiliary scales can 
then be related to the original integrals by taking limit of the auxiliary external scales da 02]. 
Unfortunately, the calculation of Feynman integral becomes more complicated as the number of 
scales increased. 

In this work, we solve single scale soft phase space integral by introducing internal auxil¬ 
iary scales. Differential equations with respect to these internal scales are derived and solved. 
Then the internal auxiliary scales are integrated over, and the desired single scale integral is 
recovered. There are several advantages in this approach. First, the internal auxiliary scales 
will be introduced by inserting a delta function. Some of the momentum integral can be in¬ 
tegrated out trivially by the delta function. Effectively, the number of loops involved in the 
calculation is reduced by one. Reducing the number of loop significantly simplifies the com¬ 
plexity of Intergration-By-Parts (IBPs) reduction [iHl 03 120] • Second, the property of scale 
invariance for eikonal integral can be fully exploited after the introduction of internal auxiliary 
scales. The consequence is that auxiliary integral only depends on a scale invariant combination 
of the auxiliary scales. Therefore the functional dependence of auxiliary integral on the auxil¬ 
iary scales has very restricted form. Third, the internal auxiliary scales are introduced in a way 
that integrating over them is also trivial. In fact, only simple interative integrals of harmonic 
polylogarithms (HPLs) [21] are needed for this calculation, in contrast to the extensive uses of 
multiple polylogarithms in Ref. [T2] . 

^We work in dimensional regularization with space-time dimension continued to d dimension throughout this 
work. 
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The single scale soft phase space integrals discussed in this work are related to the soft-virtual 
corrections for color singlet particle production at hadron collider. Perhaps the most important 
example is the soft-virtual corrections for Higgs production mm- Consider the process 

Ni{Pi) + N 2 {P 2 ) ^ H{q) + X . (2.1) 

Thanks to QCD factorization, the hadronic cross section for this process can be written schemat¬ 
ically as a convolution of parton distribution functions and partonic cross section 

(T = f ® f ® a . (2.2) 

The partonic cross section (t{z) depends on partonic threshold variable 2 ; = where Mh is 

the mass of Higgs boson, and \/5 is partonic center-of-mass energy. The soft-virtual corrections 
are referred to the leading term of a{z) in the expansion of z around 1. Since the energy of 
final state radiation is Ex ~ — z)/2 + 0{{1 — z)"^), calculation of soft-virtual corrections is 

equivalent to calculation of differential distribution in Ex, in the limit where Ex is much smaller 
than all other scales in the problem. In this limit, the leading power contributions in (1 — ; 2 ) 
come from the partonic process 


gipi) + 9{V2) ^ H{q) + X . (2.3) 

In QCD perturbation theory, X contains various combination of massless QCD parton. For 
example, at NNLO, the double-emission contribution contains the process 


gipi) + g{.V 2 ) H{q) + g{ki) + ^(fca), 
while the triple-emission contribution at N^LO contains the process 

9{Pi) + 9{P2) H{q) + g{ki) + ^(^ 2 ) + ^(fcg), 


(2.4) 


(2.5) 


where pf = kf = 0. The phase space integrals in the limit Ex Mh are single scale integrals, 
whose dependence on Ex can be determined uniquely by dimensional analysis. The main purpose 
of this work is the calculation of these integrals. 

We describe in detail the method mentioned above, taking double-emission soft phase space 
integrals as an example. After application of IBPs reduction E|, the result for double-emission 
contribution to NNLO soft-virtual corrections can be expressed in terms of two nonzero two-loop 
master integrals. 


h 

I 2 


f ^ ^+{kl) 5+{kl)5{2ExMH - {pi+P2)-{ki + k2)), 

r d'fci d-^fc2 h{kl) 5+{kl) 5{2 ExMh - (pi + P2)-(fci + k2)) 

J {ki + k2y (2 pi-ki) (2pa■ ^ 2 ) 


( 2 . 6 ) 


^It also depends on Higgs mass Mh, top-quark mass Mt, renormalization scale ga and factorization scale gp, 
which are irrelevant to our discussion. 

^We use both REDUZE 2 [22] and LiteRed [23] for IBP reduction. 
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where 7 ^ = 0.577216.... The integrals are dehned in d = 4 — 2e dimension to regulate infrared 
singularities. The dependences of Ji and I 2 on the kinematical variables Ex, Mh, pi and p 2 are 
highly constrained. From dimensional analysis, we have 

h oc {Ex MhY-^ , h oc -^{Ex MhY-^ , (2.7) 

Pl-P2 

where we have also made use of the fact that the scale transformation 


Pi-tAipi, P2-^>^2P2 (2.8) 

for the soft integrals is not anomalous, he., anomalous term of the form {prP2Y'^ is always absent. 
In the rest of this paper we will set 2ExMh = 1 for simplicity because their dependence are 
known exactly. We will also let pi-p 2 = 2, and introduce the standard notation 

Pi-k = k~^, p2-k = k~. (2.9) 

As explained above, Ji and I 2 are single scale integrals. Therefore the method of differential 
equation for Feynman integral can not be applied directly. The solution is to introduce internal 
scales by inserting an unit operator 


1 = d'‘l6'‘^(l-(k,+h)) 


( 2 . 10 ) 


into the original integrals, 

h = 


I 2 = 




TT' 


■^J^Q{1 ) S{1 - {pi + P2)-I)li{l,pi,p2) , 

^ n^2{l,Pi,P2) 

^^0(/ )5{l- {pi+p2)-l)- 


P 


where we have introduced the auxiliary integrals 

d'^fci 


Xi(/,Pi,P2) y '5+((^ - hY) , 

d'^k, 6+{kl)6+{{l-k,Y) 


l2{l,Pi,P2) 


^d/ 2-1 { 2 p^-ki){ 2 p 2 -{l-ki)) 


( 2 . 11 ) 


( 2 . 12 ) 


Integrals of the form in Eq. (I2.12p are hrst introduced in the context of fully exclusive soft 
function in ref. [23], and have been used in the calculation of threshold soft function for direct 
photon production [25] and boosted top-quark pair production [26]. We will see later that they 
are particularly convenient for evaluating soft phase space integral. 

An important feature of the auxiliary integrals in Eq. (I2.12p is that they transform in a simple 
manner under Eq. (12.Sh . The scale symmetry and dimensional analysis determine that 

MI^PuP 2 ) fi (y) , 
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(]2\{d-4)/2 

Mhpi,P 2 ) —•^ 2 ( 1 /), 


(2.13) 


where functions fi{y) depend on a single kinematical variable y = ). To determine the 

functional form of fi{y), one only needs to calculate the integrals in Eq. fl2.12p with the constraints 


/+ = 1 , 


r = 1. 


(2.14) 


The full functional dependence can easily be recovered using Eq. fl2.13p . For this reason, we will 
also denote the auxiliary integrals as 




(2.15) 


This is an important simplihcation because the number of independent scale have been signih- 
cantly reduced, but not reduced too much such that there is still nontrivial dependence on the 
scale. 

The auxiliary integrals in Eq. fl2.12p have been calculated to all order in dimensional regular¬ 
ization parameter in Ref. 


Xiix) = 


X2{x) = - 


e^7Br(l - e) 

2(1 -2e)r(l -2e) 
e"T'®r(l - e) 


X 


4er(l - 2e) 


X \l-x) " 2 ^i(-e,-e; 1 - e;x) 


(2.16) 


Results in Eq. fl2.16p can be obtained by parameterizing the phase space integral in light-cone 
coordinate and performing the phase space integral explicitly. However, due to the quadratic con¬ 
straint imposed by the delta function 5{{I — fci)^), the resulting integral is not easy to computed. 
Alternatively, they can be obtained by dispersion method [27], as along as the discontinuities 
of the integrals are well understood. While the dispersion approach is intrinsically interesting, 
we prefer not to elaborate it here and leave it for future work. Instead, we use the method of 
differential equation, explained in the next section, for the calculation of such integrals. 


2.2 Solving the auxiliary integrals by differential equation 

Method of differential equation for Feynman integral is originally proposed for loop integral only. 
It was then realized that the same method can also be applied to phase space integral with small 
modihcation The reason is that on-shell condition for phase space integral can be regarded 

as Feynman propagator for the purpose of IBP reduction or calculating derivative. 




(p2 


1 

27ri 


pz 


iO 




iO 


(2.17) 


Cut propagators are then recovered as delta function after IBP reduction. The only difference of 
phase space integral IBP reduction with conventional loop integral reduction is that integral with 
zero or positive power of should be set to zero, because they do not give rise to discontinuity. 
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The derivation of differential equation for phase space integral can then be proceeded similar to 
loop integral. Taking derivative w.r.t to x can be expressed as 


dx 2x dlf^ 

The expression on the right-hand side of Eq. fl2.18p can be reduced using IBP identities. We 
obtain 


dx X 

d^ 2 (x) (l-2e) ^ , 

dx 2x(l —x) ^ 


(1 + 

(1 -x) 


X 2 (x). 


(2.19) 


It is a straightforward exercise to solve the system order by order in e, up to undetermined inte¬ 
gration constants. The differential equation for Xi (x) is homogeneous. Its integration constants 
have to be hxed by explicit calculation. Using the result given in Eq. fl2.16p [21], Xi(x) expanded 
in e is given by 


Xi(x) 


2(1 -2e) _ 


1 — eHo^x) -|- e ( — 


3C2 


+ ) + e ( — 


7(3 


+ 7tC277o( 3^) — Tfo, 0 , 0 ( 3 ^) 


e^l - 


+ -^CsHoix) — -(2Ho^o{x) + i7o,o,o,o(a^) 


( 2 . 20 ) 


where we have only expand through to 0{e'^) for simplicity. In Eq. (I2.20p . is Riemann’s zeta 
value, H^{x) are HPLs, hrst introduced in physics in Ref. m- In our case, the required HPLs 
are dehned recursively by 


/■“dt /■“ dt 

Ho,M = / -HAt ), . ( 2 . 21 ) 

Jo ^ Jo ^ 

If w consists of O’s only, then it is dehned to be Hq^{u) = W u/n\. As a standard notation, we 
replace {k — 1) O’s followed by a 1 with k in the weight vector. For example, iLo,o,i('w) = Hoiu). 
Substituting Xi(x) into Eq. fl2.19p and solving the differential equation for X 2 (x), we obtain 


Xofx) 


-1 


4e(i _ x) + &?'77i(x)) + (hf + 6f%i(x) + i7o,o(x) - H^,o{x) 

+ ^2^- 7 fi^i(x)^ -|- ^62^ + 2 ^ 2 Ho{x) + &2 ^ 7 Li(x) -|- 63 ^- 7 fi_i(x) — Hofifi{x) + Hi^Qfi{x) 

— Hi,i,o{x) + Hi^i^i{x)^ + €^^62' -|- -( 3 Hq(x) + Hi(x) — -(2Hofi{x) + -(2 Hi^q{x) 

+ ^2^7fi^i(x) -|- b^^Hi^i i{x) + Hqq^q^o{x) — Hi^ofio{x) + Hi iqo{x) — Hi^i^i^o{x) 

+ 6™i/i,i,i,i(x)) , (2.22) 
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where b ^2 integration constants yet to be specified. An obvious way to fix these constants 
are calculating them on the boundary explicitly, x = 0 and x = 1. However, X 2 {x) is singular 
in either of these two limits, as can be seen from the known result, Eq. fl2.16p . Setting naively 
a; = 0 or X = 1 in the integrand doesn’t work and one has to perform expansion carefully along 
the line of the method of expansion by region [2S]- This however requires further study on the 
application of the method of expansion by region to phase space integral, which is beyond the 
scope of current work. Instead, we fix these constants by going to d = 6 — 2e dimension and 
requiring that X 2 (x) in 6 — 2e dimension vanishes in the limit of x ^ 0. Using dimensional 
recurrence relation [ 22 ] 0 , X 2 (x) in 6 — 2 e dimension is given by 

= ( 2 . 23 ) 

The condition hma,^o (x) = 0 gives 


b2 


2 ’ 3 ’ 16 /’ 


(2.24) 


where b ^2 ^^e {i + l)-th element of the vector 62 - We have therefore determined the nontrivial 

auxiliary integral X 2 (x) in almost an algebraic way. In particular, the integration constants are 
fixed by regularity condition alone. 


2.3 Prom auxiliary integral to soft phase space integral 

With the auxiliary integrals at hand, we are ready to solve the soft phase space integrals in 
Eq. fl2.6p and fl2.11l) . in which we are actually interested. To this end, we need to restore the full 
dependence of Xj(/,pi,p 2 ) on /, /+ and l~. 
have 

Zl{l,Pl,P2) 

l2{l,Pl,P2) 

where 

and Xi{y) and X 2 {y) are given in Eqs. (I2.20p and (I2.22p . Under light-cone decomposition, the 
original integrals in Eqs. (12. 6 p and (12.lip can be written as 

h =^^^57531 / d/- 0 (/+/- - /2) 5(1 - /+ - /-)X,(/,P1,P2). 

■^We use the implementation of these relations in the package LiteRed. 


Using dimensional analysis and scale symmetry, we 

(2pi-p2)P 


(2.27) 






















Since neither the integrand nor the delta function constrain depends on we can integrate the 
transverse angular components out explicitly 



n{d-3) 



(2.28) 


where 


27r(^+h/2 

nm) = — - 

^ ^ r((D + i)/2) 

Introducing the parameterization 

= /+r - |/^|2=M/+r, 


(2.29) 


(2.30) 


/+ and I can be easily integrated out in closed form, because Ti(y) are functions of u only, 

(2.31) 


^ 2 Pl-P 2 
^ ( 21 +) ( 2 /-) 


= u. 


For Ji the remaining integral in u can be readily carried out order by order in e. However, 
additional complexity arises for I 2 , for which the explicit integral can be written as 


h = 


/ du (/+(!-/+)) ^ - u)-^u-^X2(u). 


(2.32) 


We note that the integrand contains an unregulated singularity at u = 0. To regulate the 
integral, one approach is to extract the Inu dependence of X 2 (u) to higher order in e and resum 
it as and then apply the method of sector decomposition [30] for the integral. Alternatively, 
the singularity becomes a spurious one in 6 — 2e dimension, and we can calculate it directly 
without any regularization procedure. We can then use dimensional recurrence relation for the 
soft phase space integral to obtain the result in 4 — 2e dimension. The advantage of this approach 
is that there is no power divergence in the integrand in 6 — 2e dimension, and the integral can 
be straightforwardly carried out order by order in e. Similar strategy has been employed in the 
direct calculation of soft phase space integral [12], in the calculation of multi-box diagrams [3T] . 
and in searching for quasi hnite master integral basis [32] . In d = 6 — 2e dimension, Eq. fl2.32p is 
now given by 


h = 


e^7B(](3 _ 2e) 


TT 


2-e 


d/+ / du /+(!-/+) (1 -m) 


l-2e 


l-e.,-lT6-2e 

Xn 


u 


[U] 


(2.33) 


where we have kept the factor unchanged under dimensional shift by convention. Eq. fl2.33p 
still contains the singular factor l/u in the integrand. However, since the function X 2 ~‘^^(u) 
vanishes in the limit of u —?■ 0, the integral is hnite. We can then expand the integral order 
by order in e, and integrate in the variable u using the dehnition of HPLs in Eq. fl2.2ip and 
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integration-by-parts relations. In this way we obtain results for the original soft phase space 
integral in c? = 6 — 2e dimension, 


r6-2e 

-'1 


r6-2e 

-'2 


20160 


1 + 


293e 


+ e^ 


508636 

11025 


-rC2 +e 




27711436 _ 293C2 _ 62(3\ 


128625 


3 J 


+ e 


/113183745616 

508636(2 

18166(3 


V 121550625 

1575 

105 

4 j + P 


4254271875 


27711436(2 31535432(3 ^ 434(2(3 2637(4 


18375 33075 

;/21179542653666496 


_ 1022(5 \ 

3 140 5 J 

113183745616(2 1718109032(3 


1922(1 


9 


1 


1340095640625 
127159(4 
1225 
^ /31 

Cj + dy- 


17364375 

42778(5 2473(6 \ 

25 1 ^) 


385875 


+ 


18166(2(3 

15 


19(2 


~ 5(3 ) -|- e 


601 328(2 95(3 


9 


9 


+ 2(4 


/9451 

5308(2 

1511(3 , 

161(2(3 

, 38(4 , 4/ 

V 27 

27 

9 + 

3 

+ 3 37&j+£( 


78904(2 


81 


81 


22541(3 , 3059(2(3 ^ 295(1 ^ 2039(4 703(5 ^ 4859(6 ^ 


27 


9 


36 


48 


) 


(2.34) 


where we have kept the expansion through to transcendental weight 6. Then we can convert the 
results back to d = 4 — 2e dimension ones using the following dimension recurrence relations 


6-2e 


h 


8(-3 + 2e)(-7 + 4e)(-5 + 4e)/° 

-1 + e 

4(-3 + 4e)(-l + 4e) (-315 + 1902e - 4050e2 + 39646^ - 1824e^ + 320e^) 
' (-l + e)e3(-l + 2e) 

4(-l -F e)(-3 + 4e)(-l + 4e)/®-'^ 


(2.35) 


Substituting Eq. fl2.34p into Eq. fl2.35p . we obtain 


Ii = 


24 


1 + 


22e 


+ e" 


340 


7 ( 2 )+ £( 


4600 154(2 


27 

_ 1364(3 _ 9(4^^5^724192 _ 32200(2 


62(3 \ 
3 ) 


58576 2380(2 


9 


4 J 


243 


27 


21080(3 ^ 434(2(3 


81 

33(4 


9 

1022(5 


27 


3 


h = - 


/8822080 410032(2 

729 81 

6^(2 - 24e3(3 


285200(3 ^ 9548(2(3 ^ 


81 


9 


16e3 


23 


12 ^ V 9 


3 

1922(2 
9 

3436(5 j 


15 


' - 85(4 - 


2 5 

22484(5 


15 

6(2512( 


9 


2473(6 \ 

16 J 

12539(6y 
144 ) 
(2.36) 


These are the correct results for the NNLO soft phase space integral. We have seen in this section 
that the combined use of differential equation, integral of HPLs, and dimensional recurrence 


10 


















































































relation leads to an almost algebraic way for the determination of soft phase space integral. We 
will apply this method to the calculation of soft triple-emission phase space integrals in next 
section. 


3 Soft triple-emission phase space integral 

For soft-virtual corrections at N^LO, we need to calculate phase space integrals with three soft 
partons emitted in the hnal state. After reduction by IBP identities, we hnd a set of 7 master 
integrals need to be calculated m, 

Ji = J[dk], 

^ W+WWh+W) ’ 

^ (^1 + k2 + k^y (2/^2 ) (2(A;i -|- ^3)+) 

^ (^1 + k^y {ki + k2 + k^y i‘^kf) {2k2 ) (2(/c2 -l- 

^ {ki -|- /c2)^(/ci -l- (2(/ci -|- k2)~^) (2(/ci -|- kz)~) 

J7=j[dk] ( 2 fc+) ( 2 ^ 2 -) ’ 

where we have used the short-hand notation 

r p3e7B p 

J d%d%d%d+{kl)d+{kl) 6 +{kl)d{l - {k, + k2 + k,) • (pi +P2)). ( 3 . 2 ) 

The calculation of these integrals follows closely the method given in Sec. HI 


3.1 Auxiliary integrals for soft triple-emission phase space integrals 

As in the double-emission case, we introduce internal auxiliary scale by inserting into Eq. fl3.ip 
a delta function, 

1 = Jd^isi<i)^i_^k^ + k2 + ks)). (3.3) 

To dehne the auxiliary integrals, it is convenient to introduce the corresponding auxiliary topology 
in Fig. [U Auxiliary integral can be labeled by the power of its individual propagator, 

J{x) = Topo(ji,j2,...,j9) 
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I— 
I 



_I 

Figure 1: Auxiliary topology for the auxiliary integrals. Dash line is phase space cut, and double 
line are eikonal propagators. 


3 2e7is 


TT' 


X 


d-2 


d% d% d'^h - ih + k2 + h)) 


[k!]i^ [ki]f [kiy; 


[{k, + k2)T [{ki+ks)T mi + k2)m mi+ksrr mr mr 

We can then write the auxiliary integrals as 
Ji{x) = Topo(l, 1,1, 0, 0, 0, 0, 0, 0) 


(3.4) 


o'if-'yE 


d^ki d‘^k2 5mi)m^i)mii -ki- k2f) , 


TT' 


d-2 


J 2 {x) = Topo(l, 1,1, 0,1,1, 0, 0, 0) 


3 2e7B 


TT' 


d-2 




smmikDsm - ki - k2f) 

(/ — k2y{2{ki + ^ 2 )+) 


Jsix) = Topo(l, 1,1, 0, 0, 0,1,1, 0) 






J {2k:^){2{l - k2) + ) 


J 4 (a:) =Topo(l,l,l,0,l,l,0,l,l) 


3 2e7B 


TT' 


d-2 


d^ki d'^fca 


5mi)mki)m{i-h-k2f) 

{l-k2)\2kt){2k^){2{l-k^)+)' 


J 5 (x) =Topo(l,l,l,l,l,0,0,l,l) 
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=2e7B 


.J,, .i,. h(k^)h(k^)h((l -k,- k,f) 


TT' 


d-2 




(fci + A;2)2(/ - fc2)2(2A:2+)(2(/ - ^ 2 )") ’ 


J 6 (x)=Topo(l,l,l,l,l,l,l, 0 , 0 ) 




TT' 


d-2 




{h + A;2)2(/ - fc2)2(2(fci + fc2)+)(2(/ - ^ 2 )-) ’ 


(3.5) 


where they correspond to the original integrals Ji_ 6 - The calculation of J 7 requires additional 
care and will be addressed later in Sec. 13.21 The calculation of the auxiliary integrals in Eq. (I3.5p 
follows closely the procedure illustrated in Sec. El We present the results here and leave the 
details to appendix lAl 


Ji{x) 


J2{X) 

Mx) 


Ja{x) 


J^{x) 


X 


4(1 -2e)(2-9e-9e2) _ 
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+ 50 ( 27 ^ 0 ( 2 :) - 35(277i(a:) + 4773 (x) 
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y ^ + —CsHoix) -+ 30C2i^2(a:) - 100 ( 2 ^^ 0 , 0 ( 3 :) 
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-1 
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(3.6) 


where we have given the results through to transcendental weight 4 only. The corresponding 
expression through to weight 6 in electric form can be requested from the author. 

Our next step is to integrate these auxiliary integrals over the auxiliary scales to obtain soft 
triple-emission phase space integral. As was explained in Sec. El to further integrate the auxiliary 
integrals, we need to go do 7 = 6 — 2e dimension. Using dimension recurrence relations generated 
automatically by LiteRed, we obtain the auxiliary integrals in 6 — 2e dimension. After restoring 
the full kinematic dependence of the auxiliary integrals using scaling symmetry and dimensional 
analysis, we can straightforwardly calculate the 6 — 2e dimension soft phase space integrals using 
Eq. fl2.27p . Through to weight 6, they are given by 
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(3.7) 


Finally, using dimension recurrence relations for soft phase space integral generated by LiteRed, 
we obtain the results in 4 — 2e dimension, 
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(3.8) 


5 ’ 5 

It is interesting to note that the rational part of J 5 and Jq equals up to an overall normalization. 


3.2 The factorizable integral 

We still have the last integral J 7 to calculate. While it is straightforward to solve the correspond¬ 
ing auxiliary integral using differential equation, it is however difficult to hx the integration 
constants. In particular, the regularity condition does not hx the integration constants. To 
calculate this integral, we observe that the set of momentum {fci,/c 2 } and {^ 3 } are entangled 
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together only by the delta function, 


+ k^)-{pi + P 2 )) ■ 


(3.9) 


It is well-known in resummation study that such entanglement can be factorized into a double¬ 
emission integral and a single-emission integral via Laplace/Fourier transformation. Even more 
straightforwardly, we can consider a one parameter integral 


Jriuj) 




,du K{kl)5+{kl)5+{kl)5{u:-{ki + k2 + h)-{pi+P2)) 


d^ki d'^k2 d'^ks 


T^3d/2-3 J (fci + ^ 2)2 (2A;^) ( 2 ^ 2 ") 

The original integral is simply To factorize the integral, we rewrite it as 

= / dcui </?,12(1^1) / dci;2 </?,3(1^2) dioj — oji — UJ2) y 


(3.10) 


(3.11) 


where 
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3 2e7B 
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d^ki d^k2 
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77' 
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d'^k:iS+(kl)S{(j - /(3 ■ (pi +P 2 )) 


(3.12) 


We note that by dimensional analysis, ^ 7 , 12 ( 1 ^) is proportional to I 2 in Eq. fl2.36p . In d = 4 — 2e 
dimension, we have 


.-l-4e T 
12 


d7,12(<^) — ^ J-2 

and 77 , 3 ( 0 ;) is the simple soft single-emission phase space integral. 


(3.13) 


77 , 3 ( 1 ^) — 
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l-2e 


■ 3 2 ^ 7^3 15^4 7 ( 2(3 31C5\ , ,6/49(1 399C6\ 

2(1-2e) L 2 + " 1^8 ) 


(3.14) 


Substituting Eqs. (I3.13p and (I3.14p into Eq. (13.lip , we obtain 
1 - — + e" 
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V 10 


2(2\ , 3/ 126 

, 77(2 

10(3 ^ 
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) 

> 4(2(3 1 

2387(4 
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71C3 3 + 

30 

3 / 

+ e° 
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36(2 


+ 84(3 


154 ( 2(3 10(1 2201(4 , 2849(5 1111(6 


15 


10 


+ 


15 


18 


(3.15) 


Eqs. (13.8p and (I3.15P are the main results of this work. Similar integrals have also been 
computed in Ref. jT2] using completely different method and our results agree perfectly with 
theirs after taking into account the difference in the normalization. 
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4 Conclusion and outlook 


In this work we have presented the calculation of soft triple-emission phase space integrals, which 
are relevant for soft-virtual corrections to Higgs production at N^LO [TOl ITT] . The major tech¬ 
niques used in this calculation are differential equation for auxiliary integrals, and integration 
of harmonic polylogarithms. The essential idea is to introduce internal scales into single scale 
integrals. Since integral in the eikonal limit enjoys the property of scale invariance, the internal 
scales should be introduced in a way that doesn’t spoil the scale invariance. In this way the 
resulting auxiliary integrals obey strong constraints from scale invariance, and therefore simplify 
signihcantly the solution of differential equations. A particular nice feature of the system of 
differential equation is that, all the nontrivial integration constants can be determined by regu¬ 
larity condition in 6 — 2e dimension. Furthermore, the solution of differential equations, as well 
as the subsequent integration of auxiliary scales, only involves the well-known generalization of 
logarithms and polylogarithms, the harmonic polylogarithms. Thanks to all these techniques, 
we are able to calculate the soft phase space integrals in an almost algebraic way. The meaning 
of such a calculation is two folds. First, it provides an alternative method to the one used in the 
calculation of Ref. [12]. Second, it provides the hrst independent conhrmation of the integrals 
calculated in Ref. [T2] . 

Recently, signihcantly progresses have been made towards the calculation of N^LO QCD 
corrections to the Higgs boson production via gluon fusion [SS] [Ml [35l [Ml EH jSH |39l EOl IHl W2\ 
iiHiliiniiiaiigiiHiiHiiiniiiaEoiisiiEHisa]. We hope that the method and results presented 
in this work will be useful in the pursuit of ultimate precision for Higgs physics. 
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A Solving difFrential equations for auxiliary integral 

In this appendix we solve the auxiliary integrals in Eq. fl3.5p using the method of differential 
equation. It turns out that the system of auxiliary integrals in Eq. fl3.5p is not closed under the 
Additional auxiliary integrals are needed. A complete master integral basis is 

(2e-l)(3e-2)(3e-l)Ji(x) 

X 

- e^(2e - l)J^{x ), 

- e^{x - l)Ji{x ), 


derivative of x. 
given by 

9i{x) = 
92{x) = 
9^{x) = 
9^{x) = 
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(2e-l)(36 - 2)(36 -l)Ji(x) 
4x 


g^{x) = - e^{x - l)xj^{x ), 
gQ{x) = - e^(x - l)xJe(x), 


grix) = e^j 2 (x) - 2e^Topo(l, 1,1, 0,1,1, 0, 0,1) - 


98(x) 
99 (x) 
9io(x) 
9ii(x) 
912 (x) 

9i3(x) 


e3Topo(l,l,l,0,0,0,0,l,l), 
e^Topo(l, 2,1, 0, 0,1,1, 0, 0), 
6='Topo(l,l,l,0,0,l,l,0,0), 

- e^(x - l)Topo(2,1,1, 0, 0,1,1, 0, 0), 
e^xTopo(l, 1,1,1,1,1,1, —1, 0), 


_ 6e^Topo(l,l,l,0,0,0,0,1,1) 3 

X — 1 

^ 3(26 - 1)(36 - 2)(36 - l)Ji(x) 

4(x — l)x 


xTopo(l, 1,1,1,1, 0, 


- 1 , 1 , 1 ) 


(A.l) 


Note that we have chosen a basis such that the system of differential equation has optimal 
form [SI IS]- The advantage of the optimal basis is that the resulting differential equation can 
be solved almost trivially, up to some integration constants, as demonstrated in many non-trivial 
examples [SI ETJ |S1 EHl [131 EQl EB EB uni ISl ESI ISl EB EB ESI EB E2] • Specifically, in d = 4 — 2 e 
dimension, the differential equation for the optimal basis can be written as 


dgjx) 

dx 


6 



1 



9{x). 


(A. 2 ) 


The important feature of optimal basis is that the dependence on e can be factored out, and 
A and B are constant matrices in field of rational numbers. Since the e dependence on the 
right-hand side of Eq. (IA.2P has been factored out, it is clear that the derivative of g{x) at 
depends only on the 0(e*) terms of g{x). Therefore, one can simply solve Eq. flA.2p in a 
bottom-up way by integrating the right-hand side of Eq. (IA.2D . In our case, the matrices A and 
B are given by 
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5 
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-I 0 0 
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(A.3) 
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and 


B 


/ 0 0 0 0 0 
0 0 0 0 0 
-i 0 2 00 
I 0-320 
I 0 0 0 2 
0 0 0 0 0 
0 2 6 0 0 
i 0 0 0 0 
-^0 0 00 
0 0 00 
0 0 0 0 0 
-^0 0 00 
VI 0 0 0 2 


0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 6 0 0 

2 0 0 0 0 

0-10 00 

0 0 2 0 0 

0 0 0 -| f 

0 0 0 -| I 

0 0 0 0 0 

2 0 0 -| I 

0 0 18 0 0 


0 0 0 \ 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

1 0 0 

0 0 0 

0 0 0 

-loo 
- 0 0 

3 0 0 

-I -1 0 

0 0 -1 / 


(A.4) 


While integrating the differential equation is trivial with the optimal basis, determining the 
integration constants is not. In particular, the master integrals at phase space point x = 0 or 
X = 1 are in general singular in 4 — 2e dimension, therefore method of asymptotic expansion is 
needed. However, for phase space integral, we are not aware of simple procedure for asymptotic 
expansion. The solution of this difficulty is by going to d = 6 —2e, where the degree of singularity 
is lessen. In particular, we hnd that integration constants for all the master integrals, except the 
ones with trivial x dependence, J7i(x) and 1 / 2 (x), can all be uniquely determined by the following 
regularity conditions. 


• There are no power law divergences for the master integrals at x = 0 nor x = 1. 

• If the master integrals in the x 0 limit are free of logarithmic singularity, then they 
vanish in that limit. 


These conditions can be inferred by studying the asymptotic expansion of corresponding auxiliary 
integrals, treating the cut propagators as normal propagators. For the integration constants of 
J7i(x) and J72(x), we calculate directly using the dispersion method, to which we wish to come 
back in the future. 
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